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ABSTRACT 



lO Spectra of exoplanet atmospheres provide us the opportunity to improve our understanding of 

,_k these objects just as remote sensing in our own solar system has increased our understanding of the 

(^— ^ solar system bodies. The challenge is to quantitatively determine the range of temperatures and 

ry\ species abundances allowed by the data. This challenge is often difficult given the low information 

y-^ content of most exoplanet spectra which commonly leads to degeneracies in the interpretation. A 

H^' variety of temperature and abundance retrieval approaches have been applied to exoplanet spectra, 

but no previous investigations have sought to compare these approaches. In this investigation we 
compare three different retrieval methods: Optimal Estimation, Differential Evolution Markov Chain 
Monte Carlo, and Bootstrap Monte Carlo. We call our suite of retrieval algorithms the CaltecH 
C^ Inverse ModEling and Retrieval Algorithms (CHIMERA). We discuss what we can expect in terms 

of uncertainties in abundances and temperatures given current observations as well as potential future 
observations and what conclusions can be drawn given those uncertainties. In general we find that the 
three approaches agree for high quality spectra expected to come from potential future spaceborne 
missions, but disagree for low quality spectra representative of current observations. We also show that 
the Gaussian posterior probability distribution assumption made in the Optimal Estimation approach 
is valid for high quality spectral data. Furthermore we compare the results from a parameterized 
temperature profile versus a full classical Level-by-Level approach and discriminate in which situations 
each of these approaches is applicable. We also discuss the implications of our models for the inferred 
C to O ratios of exoplanetary atmospheres, which of course are important for understanding formation 
environments. More specifically we show that in the observational limit of a few photometric points, 
the retrieved C/0 is biased towards values near solar and near one simply due to the assumption of 
uninformative priors. 
Subject headings: planetary systems — planets and satellites: atmospheres — radiative transfer — 
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1. INTRODUCTION 

Thermal em.ission spectra (^ 1-30 microns) of extraso- 
lar planets can tell us about their atmospheric temper- 
atures and compositions (see e.g., Tinetti et al. 2007; 
2010; Grillmair et al. 2007; 2008; Swain et al. 2009a; 
2009b; Madahusudhan & Seager 2009; Madhusudhan et 
al. 2011; Lee et al. 2012; Line et al. 2012). At the 
moment, these observations come in two types, broad- 
band photometry mainly from the Spitzer Space Tele- 
scope (see e.g., Knutson et al. 2010) and ground-based 
instruments (CroU et al. 2010; Anderson et al. 2010; 
Gibson et al. 2010; Deming et al. 2012; Gillon et al. 
2012 ), as well as higher resolution spectra such as Hubble 
Space Telescopes Wide Field Camera 3 (WFC3) (Berta 
et al. 2012; Swain et al. 2012; Deming et al. 2013 ) 
and Near Infrared Camera and Multi-Object Spectrom- 
eter (NICMOS) (Swain et al. 2009a; 2009b; Tinetti et 
al. 2010; Gibson et al. 2011; Crouzet et al. 2012). From 
these observations, signatures of a variety of molecules 
have been detected including H2O, CH4, CO and CO2 
(Tinetti et al. 2007; Swain et al. 2009a; 2009b; Tinetti 
et al. 2010), although the robustness of some of these de- 
tections have recently been called into question (Gibson 
et al. 2011). These same data have been used to infer 
the presence of atmospheric temperature inversions for a 
subset of hot Jupiters (e.g.. Burrows et al. 2007, Knut- 
son et al. 2008;2010; Forntey et al. 2008; Madhusudhan 
& Seager 2009; 2010). 

While the above studies have given us insight into the 
nature of these planetary atmospheres, very few have fo- 
cused on the uncertainties in temperatures and compo- 
sitions. Until relatively recently (Madhusudhan & Sea- 
ger 2009, Madhusudhan et al. 2011, Lee et al. 2012, 
and Line et al. 2012), most compositions and temper- 
atures and thus the subsequent conclusions, were de- 
termined through self-consistent forward modeling ap- 
proaches that only explore a few potential solutions with- 
out a well-defined characterization of the uncertainty dis- 
tributions of the physical parameters (e.g., Burrows et al. 
2005; Fortney et al. 2005; Burrows et al. 2007). Fur- 
thermore, some self-consistent solutions make physical 
assumptions that may not necessarily be valid in exo- 
planetary atmospheres such as the assumption of ther- 
mochemical equilibrium gas concentrations or radiative- 
convective temperature structures (that is they may ig- 
nore other potentially important processes such as verti- 
cal mixing, photochemistry, zonal winds etc.). Addition- 
ally, this forward modeling approach is often unguided 
by the data and primarily driven by preconceived no- 
tions of how the atmosphere "should" look (as pointed 
out by Lee et al. 2012 and Benneke & Seager 2012) with 
the best solutions being the few that provide the lowest 
values of chi-squared. 

In order to more rigorously characterize the ranges 
of allowable temperatures and compositions, Madhusud- 
han & Seager (2009) developed a multidimensional grid 
search approach which can fully characterize the un- 
certainty distributions for each parameter. Subsequent 
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studies (Madhusudhan et al. 2011; Benneke & Seager 
2012) used the more sophisticated Markov Chain Monte 
Carlo approach (MCMC) to accomplish this goal. How- 
ever, such approaches require the computation of many 
millions of models in order to fully characterize the pa- 
rameter uncertainties which may not be feasible for more 
sophisticated forward models with many free parameters. 
In order to remedy this problem Lee et al. (2012) and 
Line et al. (2012) used the much faster Optimal Estima- 
tion (e.g, Rodgers 2000) approach to estimate the error 
distributions of each parameter. This approach is much 
faster due to the assumption that the parameter error 
distributions are Gaussian. However, this Gaussian as- 
sumption may result in an incorrect estimate of the error 
distributions (Benneke & Seager 2012). 

The goals of this paper are to first understand the com- 
position and temperature uncertainty distributions for 
different degrees of observational quality, and second to 
understand how those derived uncertainty distributions 
differ between the two fundamental parameter estima- 
tion approaches. Optimal Estimation and MCMC. This 
investigation represents the first direct comparison and 
synthesis of these retrieval approaches as applied to ex- 
oplanet atmospheres. A secondary goal is to understand 
how the derived composition uncertainties propagate into 
the C/0 uncertainty. We accomplish these goals by com- 
paring three different retrieval algorithms: Optimal Esti- 
mation (OE), a new MCMC algorithm known as Differ- 
ential Evolution Markov Chain Monte Carlo (DEMC), 
and the model-dependent Bootstrap Monte-Carlo ap- 
proach (BMC). This investigation is analogous to the 
investigation carried out by Ford (2005) on radial ve- 
locity data. First we will describe the three different re- 
trieval techniques as well as our forward model in fp] We 
call our three-pronged retrieval approach CHIMERA- 
CaltecH Inverse ModEling and Retrieval Algorithms. 
Second, we compare the three spectral retrieval meth- 
ods on different synthetic spectral data sets of varying 
observational quality in order to asses the robustness of 
the error estimations from each approach in fpj We will 
also compare the parameterized temperature profile ap- 
proach (e.g., Madhusudhan & Seager 2009; Line et al. 
2012) with the Level-by-Level profile approach (Lee et 
al. 2012). Finally, we will discuss the implications of 
these uncertainties for the estimated C to O ratios. 

2. METHODS 

In this section we describe the retrieval techniques, the 
forward model, and the parameterizations we use to re- 
trieve the temperatures and compositions from thermal 
emission spectra. 

2.1. The Retrieval Techniques 

We use three different retrieval techniques to infer the 
compositions and temperatures from a spectrum. The 
techniques are inherently Bayesian as they attempt to 
solve the inverse problem by summarizing the full shape 
of the posterior in terms of the location in parameter 
space of the maximum likelihood and the uncertainties 
about that location. The first, and the fastest (least 
number of forward model calls) of these approaches, is 
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Optimal Estimation, the second is the model-dependent 
Bootstrap Monte Carlo, and the third is Differential Evo- 
lution Markov Chain Monte Carlo. 

2.1.1. Optimal Estimatton (OE) 

The Optimal Estimation retrieval approach is well es- 
tablished in the fields of Earth atmosphere remote sens- 
ing (Rodgers 1976; Towmey 1996; Rodgers 2000; Livesay 
et al. 2006; Kuai et al. 2013), solar system atmosphere 
remote sensing (Conrath et al. 1998; Irwin et al. 2008; 
Nixon et al. 2007; Fletcher et al. 2007; Greathouse et al. 
2011), and recently exoplanet atmosphere remote sensing 
(Lee et al. 2012; Line et al. 2012). The basic approach 
is to minimize a cost function to obtain the maximum 
a posteriori (MAP) solution. Using Bayes theorem and 
the assumption that the data likelihood, the prior, and 
the posterior are Gaussian , one can derive the following 
cost function (or log likelihood): 



J(x) = (y-F(x))^S-i(y-F(x)) 

+ (x-Xa)^S-l(x-Xa) 



(1) 



where y is the set of n observations, x is the set of m 
parameters which we wish to retrieve or the state vector, 
F(x) is the forward model that maps the state vector 
onto the observations, and Sg is the n x n data error 
covariance matrix (typically off diagonal terms are zero 
and the diagonal elements are the square of the la er- 
rors of the observations). Xa is the a priori state vector 
and Sa is the m x m a priori covariance matrix. The 
first term in equation [T] is simply "chi-squared" and the 
second term represents the prior knowledge of the param- 
eter distribution before we make the observations. For 
high quality observations the second term is generally 
not important as most of the information in constraining 
the state vector comes from the observations. For low 
quality observations it is just the opposite. Following Ir- 
win et al. 2008 we minimize equation [T] with Newton's 
iteration method given by 

xi+i = Xa + Sa'K?'(KiS-iK?' + S-1) 

(F(x)-y-Ki(xa-Xi)) (2) 



where i is the iteration index and Kj is the Jacobian 
matrix at i {Knm — §^)- Rather than taking the full 
Newton step we damp the solution with 



^i+i 



Xi 



Xj+l - Xj 

1 + C 



(3) 



where C is the damping parameter. At each iteration we 
evaluate J(xi+i) and J(x;_|_]^). If the latter is smaller, 
we set the state vector for the next iteration to Xj , -^ and 
decrease ^ by 0.3. Otherwise we keep increasing C by a 
factor of 10 and re-evaluate equation [3] and equation IT] 
until J(x;^-^) becomes less than J(xi+i). Convergence is 
achieved when J changes by less than 1 x 10^® from the 
previous iteration, which typically occurs after '-^lOs of 
iterations. The resulting state vector is the MAP solu- 
tion. The uncertainties on the state vector parameters 
are given by the posterior covariance matrix 



S = (KTSg^K + Sa^)-^ 



(4) 



This matrix represents a multi-dimensional normal dis- 
tribution (see Rodgers 2000 for the derivation). The di- 
agonal elements are the square of the marginalized er- 
rors whereas the off diagonal terms describe the corre- 
lations/degenracies amongst the parameters. The first 
term ,K'^Sg """K, represents the uncertainties due to the 
measurement errors. This term uses the local gradi- 
ent information to estimate the parameter uncertainties. 
The second term represents the prior uncertainties be- 
fore making the measurements, which has less influence 
for higher quality data. Again, the major assumption in 
equations [l] and [4] is that the posterior for each paramter 
is Gaussian. This assumption is only valid when the re- 
gion in phase space over which the forward model can be 
linearized is broader than the parameter uncertainties. 
However, it is this assumption that allows this approach 
to be extremely fast requiring only tens of forward model 
calls. As we shall see in ^ this assumption is valid 
for data that is of "good" quality, but breaks down for 
"poor" quality data. 

2.1.2. Model Dependent Bootstrap Monte-Carlo (BMC) 

A common way to more robustly characterize errors is 
through a Montc-Carlo resampling of the data (see e.g.. 
Press et al. 2002 Chapter 15.6, Ford 2005 §4.2) in which 
many thousands of realizations of the original data (in 
our case, the spectra) are created using the uncertainties 
from the original dataset. These synthetic data are then 
refit using, say, OE, and the resulting best fit parameter 
distributions represent the uncertainties. There are mul- 
tiple ways of generating the synthetic data realizations. 
The most common way is the residual resampling ap- 
proach in which data realizations are created by adding 
the randomly shuffled residual between the best fit model 
and the data back to the original best fit model. This 
new realization is then fit and the process is repeated 
many times. The approach we take is similar, but rather 
than generate a new spectrum using the residual, we sim- 
ply take the best fit, from OE, and then resample each 
point by drawing it from a normal distribution with a 
mean given by the best fit value and the width given 
by the data error bar for that point. We chose this ap- 
proach over the residual resampling approach because 
sparse coverage spectra, like those from broadband ob- 
servations, have virtually no residual as they can be fit 
perfectly due to the greater number of parameters than 
data points. We typically generate ~1000 spectra re- 
alizations that are then refit by optimal estimation to 
obtain the state vector parameter distributions. 

2.1.3. Differential Evolution Markov Chain Monte Carlo 

(DEMC) 

The MCMC approach has revolutionized parameter es- 
timation and error analysis in many fields. It is routinely 
used in radial velocity (Ford 2005 ) and transit light curve 
(e.g., Eastman, Gaudi, & Agol 2013) error analysis. Re- 
sults from a well converged MCMC analysis can gener- 
ally be considered as the best possible representation of 
the parameter uncertainties. Recently, this approach has 
been applied to the exoplanet atmosphere retrieval prob- 
lem (Madhusudhan et al. 2011; Benneke & Seager 2012). 
Unlike optimal estimation, MCMC approaches make no 
assumptions about the shape of the posterior, but rather 
evaluate the posterior with millions of samples. 
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The basic approach of MCMC is to sample the poste- 
rior through a random walk process. The random walk 
is carried out by drawing states from some proposal dis- 
tribution and evaluating whether or not the proposed 
state has an increased likelihood over the previous. Typ- 
ically the proposal distribution is a normal distribution 
with a mean given by the current state in the chain (xi) 
and a user defined width to achieve a particular accep- 
tance rate (Gibbs sampling or Metropolis-Hastings). If 
the proposed state (xp) has an improved likelihood over 
the current state, then that state is kept (xj+i = Xp) 
and a new proposal is made from that location. If the 
proposal state has not improved the likelihood, then that 
state is either rejected or accepted with some probability. 
This previous state dependent random walk constitutes a 
Markov Chain. Given enough samples this Makov Chain 
will converge to the target posterior (see Ford 2005 for a 
more detailed explanation). 

Rather than standard MCMC approaches, we use 
an adaptive algorithm known as Differential Evolu- 
tion Markov Chain Monte Carlo (DEMC) (ter Braak 
2006; ter Braak & Vrugt 2008). The purpose of this 
approach is to obtain more appropriate proposal states 
by identifying the proper scale and orientation of 
the current estimate of the posterior. This scale and 
orientation information comes from the chain history. 
This approach gives a more efficient probing method for 
highly correlated parameter spaces and yields improved 
convergence rates. Our DEMC procedure is as follows: 

1. Apply the OE technique to the measurements to 
obtain the best fit state vector and posterior covariance 
matrix, S. This step provides an initial estimate of the 
posterior. 

2. Initialize N^nit links (xi=o-Ni„it) in each oi N chains 
(typically 3 chains, more chains will slow convergence) 
independent chains (arrays) by randomly drawing state 
vectors from the multivariate normal described by the 
posterior covariance matrix from step 1. Set the last 
link in one of the chains to the best fit state vector 
obtained in step 1. This step provides a good starting 
history from which our initial proposal states can be 
drawn. Combine each of the independent chains into 
one long chain that composes the history, Xhistory 

3. Evaluate the cost function, J, in equation [l] for the 
last link in each of the chains. If using a flat prior ignore 
the second term. This, again, is simply the equivalent of 
evaluating chi-squared. 

4. Draw two random numbers, Ri and R2, between 
zero and Nchams x *, where i is the current state in the 
chain. Initially, i = Nmu- Evaluate the proposed jump 
state given by 



Xp = Xi + 7(xRi - XR2) 



(5) 



where xr,^ and xrj are the states from different points 
in the the chain history, Xhistory 7 is a scale factor 
typically set to 2.38/-\/(2 * m) (ter Braak 2006), where 
m is the number of parameters. This factor is meant 
to give acceptance rates of ~0.23 for large m. e is a 
vector drawn from a multivariate normal distribution 



with a small variance relative to the chain variance in 
order to introduce a small amount of additional random- 
ness. Repeat this process for the other N chains — ^ chains. 

5. Evaluate the Metropolis (Metropolis et al. 1953) 
ratio, r = P(xp)/P(xi) = e~2(-'(''p)~J(^i)). If 7- is larger 
than 1, set Xi^i = Xp and if it is smaller only accept 
if it is larger than a random number between 1 and 
0. Otherwise do not update the chain, set Xi+i = Xi. 
Repeat for the other N chains — 1 chains. Add the 
updated links in all N chains to Xhistory 

6. Repeat steps 4 and 5 until convergence is met. Con- 
vergence can be determined by looking at the trace plots 
of Xhistory for cach parameter or by using the Gelman- 
Rubin statistic on the set of Nchains chains. For this we 
use the algorithm from Eastman, Gaudi, & Agol 2013 
which requires the Gelman- Rubin statistic to be less than 
1.01 and the number of independent draws to be greater 
than 1000 for each parameter . Convergence typically oc- 
curs in less than 10^ links in each of the Nchains for a total 
of Nchains X 10^ liuks. This is about an order of mag- 
nitude less than parallel tempering or pure Metropolis- 
Hastings. 

2.2. The Forward Model 

The forward model, F(x), is the most important part 
of any retrieval algorithm. It is what maps the state 
vector of retrievable parameters onto the observations. 
In the case of atmospheric retrieval, the forward model 
takes temperatures and compositions and generates a 
model spectrum. Our particular forward model numeri- 
cally solves the planet-parallel thermal infrared radiation 
problem for an absorbing, emitting atmosphere (we ne- 
glect scattering). We first divide the atmosphere into 
Nicy discretized atmospheric layers. The absorption op- 
tical depth for the fc*'' gas in the z*'' layer at wavelength 
A is 

^Tk^z,\ = fk.z^k^z.X — (6) 

IJ^atrnd 

where fk,z is the volume mixing ratio of the k*^ gas in 
the z*'' layer , (Jk.z,\ is the absorption cross section per 
molecule of the fc"' gas in z*'* layer at wavelength A, AP^ 
is pressure thickness of the z*'' slab, ^iatm is the mean 
molecular weight of the atmosphere and g is the gravity. 
The absorption cross sections are pre-computed on a 1 
cm~^ wavenumber grid at 20 evenly spaced temperature 
and log-pressure points from 500-3000K and 50-10~^ bars 
respectively (similar to Sharp & Burrows 2007). The 
cross sections for each wavelength on the pre-computed 
grid are interpolated to the atmospheric temperatures 
and pressures in the z*'* slab. To compute the total slab 
optical depth we sum the contribution from each gas to 
obtain 

Ar^,A = Y^ ATk^zA (7) 

fc=i 

Upon computing the optical depths at each level, we can 
now solve for the upwelling irradiance with 



/a = 



E 

z=0 



i?A(T,)e-^"- ^^-^Ar,„ 



(8) 
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where Ni^y is the num.ber of atm.ospheric levels and 
B\(Tz) is the Planck function at wavelength A and tem- 
perature in the z*'' slab. We use 90 atmospheric layers 
to compute the upwelling flux. 

An important part of the forward model when using 
the Optimal Estimation approach is the computation of 
the Jacobian, or the sensitivity to the upwelling irradi- 
ance with respect to the desired retrievable parameters. 
When possible, it is preferable that the Jacobian be cal- 
culated analytically for both improvements in speed and 
in accuracy. We are interested in the retrieval of both 
abundances and temperatures so we must compute the 
Jacobian with respect to both the abundances and tem- 
peratures. We make the assumption of vertically uniform 
gas mixing ratios throughout the atmosphere and hence, 
fk^z is independent of z. We now differentiate equation 
[8] with respect to the uniform gas mixing ratios for each 
gas fk to obtain 



Ojk —f, Jk 



z,\ 



2=0 



Ni^ 



Ni^ 



J2 {Bx{Tz)e~^^-^ ^"-^Ar.a ^ 



Arfc. 



J:A 



z=0 



2=z 



fk 



) (9) 



The first term is due to the changing emissivity of the 
emitting slab and the second term is how the change in 
transmittance affects the upwelling irradiance. 

The sensitivity of the irradiance to a change in tem- 
perature in the z*'' slab is given by 



dh 



(e 






-E~lrAr„. 



dB^jTz 
' dTz 



(10) 



This equation is similar to equation 14 in Irwin et al. 
2008 but we have neglected the first and last terms in 
their formula as they are small. 

Since the observations are reported as the ratio of the 
planet flux to the stellar flux and not the irradiance, we 
perform a disk integration of equations [8l |10| using four 
point Guassian quadrature and then divide by an inter- 
polated PHOENIX steUar flux grid model (Allard et al. 
2000) . 

We include only CH4, CO2, CO, H2O, H2, and He in 
our model. H2 and He are fixed in our models at ther- 
mochemically justifiable abundances. The exact abun- 
dances of these species is not critical as the sensitivity 
of the spectrum to H2 and He is minimal. We retrieve 
only CH4, CO2, CO, and H2O. We choose these species 
because they are the most spectroscopically active and 
abundant species. Admittedly we could/should include 
every possible atmospheric constituent but this would 
be unwieldy and reliable high temperature absorption 
line lists only exist for a few. On that note, we use 
the HITEMP database (Rothman et al. 2010) to com- 
pute the tabulated cross-sections for CO2, CO, and H2O 
and the STDS database for CH4 (Wenger & Champion 
1998). Below 1.7 ^m for CH4 we simply use the HI- 
TRAN (Rothman et al. 2009) database for lack of any- 
thing better (to the best of our knowledge). We use the 
Barysow et al. (2001;2002) and J0rgensen et al. (2000) 
databases for the computation of the H2-H2/He collision 
induced opacities. The Reference Forward Model (RFM- 
http://www.atm.ox.ac.uk/RFM/ ) was used to compute 



the tabulated cross sections from the line strength 
databases. We have validated our forward model through 
a detailed comparison with the Oxford NEMESIS group 
(e.g., Lee et al. 2012) and our results agree to better than 
5% (see Figure Il| (see online supplementary material). 

An additionarcomponent of the forward model is the 
instrumental function used to convolve with the high- 
resolution model spectrum. For the broadband points 
we simply integrate the flux from the high-resolution 
model spectrum with the appropriate filter function for 
that point. When fitting higher resolution observations, 
the instrumental function is assumed to be a Gaussian 
(valid for grating spectrometers) in wavelength with a 
full width at half maximum (FWHM) determined by ob- 
servations. 

Now that we have a well defined forward model we can 
define our state vector. Again, we wish to retrieve the 
abundances of CH4 , CO2, CO, and H2O and the tem- 
perature profile. More specifically, we choose to retrieve 
the log of the abundances as they can vary by orders of 
magnitude and to prevent negative mixing ratios. Our 
state vector is given by 

X = [log{fH2o)Jog{fcHj,log{fco)Jog{fco2),Tf- 

where the fkS are all assumed constant with altitude. 
We feel this is appropriate for two reasons. First, vertical 
mixing will smooth out the mixing ratio profiles over the 
thermal infrared photosphere (Line et al. 2010, Moses 
et al. 2011, Line et al. 2011), and secondly current ob- 
servations simply do not have the information content to 
warrant the retrieval of vertical mixing ratio information 
(see Lee et al. 2012). In the next section we describe 
how to go about retrieving the temperature profile. 

2.2.1. Parameterized vs. Level-by-Level (Level-hy-Level) 
Temperature Profile 

We employ two approaches to retrieve the tempera- 
ture profiles. The first, and the most commonly used in 
Earth and solar system atmosphere retrieval problems, is 
the Level-by-Level approach. This is the approach used 
in Lee et al. 2012. The second is a parameterized tem- 
perature profile approach similar to the approach used in 
Madhusudhan & Seager 2009 and Line et al. 2012. Each 
has its advantages and disadvantages described below. 

The Level-by-Level temperature retrieval approach 
seeks an estimate of the temperature at each of the Ni^v 
model layers. This approach is advantageous in that 
there are no pre-conceived assumptions made about how 
the atmospheric temperature should be parameterized. 
If the spectral signal-to-noise and resolution are high 
enough, there is generally enough sensitivity to obtain 
information about the temperature in individual atmo- 
spheric layers. However, there is a finite vertical resolu- 
tion given the quality of the observations. Typically this 
resolution is set by the width of the thermal emission 
weighting functions and how much they overlap. Gener- 
ally, when the spectra are noisy the Level-by-Level ap- 
proach fits the noise which results in unphysical structure 
in the retrieved temperature profile. This is analogous 
to fitting a high-degree polynomial to only a few points. 
There are ways to smooth unphysical structure, one of 
them being to assume a correlation among the atmo- 
spheric layers (Rodgers 2000, Irwin 2008) implemented 
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through the prior covariance matrix, Sa, with 



^a,ij y^a.ii^ajj) ' 



-ll,^(Pi/Pi)l 



(12) 



Here Pi and Pj are the pressures at the i*'' and :/*'* levels, 
respectively and h is the correlation length that controls 
the level of smoothing. The correlation length can be 
thought of as the number of scale heights over which 
the temperatures are correlated. For our simulations we 
choose h = 7 as this provides a sufficient level of detail 
without producing unphysical oscillations. When using 
this app roach our state vector is exactly as it is in equa- 
tion [TT] with T being an Ni^v vector of temperatures at 
each level. The Level-by-Level approach is only appro- 
priate when the information content of the spectra is 
sufficiently high such that the addition of the Ni^v addi- 
tional parameters is justified. For most current exoplanet 
spectra, this is an invalid approach. 

The second temperature profile retrieval approach 
makes use of a parameterization. This approach is ad- 
vantageous when the information content of a spectrum 
is low as the number of free variables is greatly reduced. 
However, the parameterization does force the retrieved 
atmospheric temperature structure to conform only to 
the profile shapes and physical approximations allowed 
by that paramterization. For our particular parameter- 
ization, we assume the atmosphere to be in radiative 
equilibrium based upon the analytic radiative equilib- 
rium temperature profile of Guillot 2010 (and others such 
as Hansen 2008, Robinson & Catling 2012). This is the 
same parameterization used in Line et al. 2012. This pro- 
file assumes two independent downwelling visible chan- 
nels of radiation and one upwelling stream of thermal 
emission. Briefly, the temperature as a function of the 
thermal optical depth, t, is given by 
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where 71 and 72 are the ratios of the Planck mean opac- 
ities in the visible streams to the thermal stream and 
the parameter a (range to 1) partitions the flux be- 
tween the two visible streams. £^2(7''') is the second or- 
der exponential integral function. The internal heat flux 
is parameterized by the temperature ,Tj„( , which is flxed 
since it has little impact on the spectra. The stellar in- 
put at the top of the atmosphere is represented by Tirr 
given by 

T,„ = l3{^y/'T^ (15) 

where i?* and T^, are the stellar radius and temperature, 
a is the semi-major axis, and /3 is a catch all term on the 
order of unity for the albedo, emissivity, and day-night 
redistribution. The grey infrared optical depth can be 
mapped onto pressure coordinates using 
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and Kjn the Planck mean thermal infrared opacity. This 
T — P mapping assumes a linear relation between the op- 
tical depth and pressure, or a pressure independent in- 
frared opacity. More complicated mappings that account 
for the pressure dependence of nm can also be used (see 
e.g., Robinson & Catling 2012). 

This temperature parameterization has 5 free parame- 
ters governing its structure: Km, 7^^, 71,3, /3, and a. Our 
parameterized state vector again, is given by equation 11 
but with T replaced with [k/a:, 7«j, 71,5, /?, a]. Com- 
bined with the gases this gives a total of 9 free param- 
eters. The temperature profiles are then reconstructed 
from the probability distributions of those 5 parameters. 

We should note that currently most exoplanet spec- 
tra often have fewer measurements than desired state 
variables. This means that each parameter cannot be 
uniquely determined. This is not a new problem (e.g., 
Madhusudhan & Seager 2009). This is why the prior is 
crucial. We can think of the prior as an "artificial" set of 
data on which the retrieval (all retrieval approaches) can 
rely when the measurements are insufficient to constrain 
a given parameter. Therefore the resulting constraints 
on a given parameter are a combination of the informa- 
tion obtained from the spectra and the prior knowledge. 
In the extreme case of no observational constraint, the 
posterior will simply be the prior. Hence, it is critical 
to choose an appropriate prior, especially for cases when 
there are more parameters than measurments. 

With the optimal estimation formalism, we can assess 
the degree to which the constraint comes from the mea- 
sured spectra versus the prior through what is called the 
averaging kernel. The averaging kernel is an rn, x tti ma- 
trix whose elements are given by 
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j^true 



(17) 



where P is the pressure, g the surface gravity (at 1 bar). 



where Xi^retrieved IS the retrieved value of the r parame- 
ter and Xj^true is the true value of the j'*'' parameter. The 
diagonal elements tell us how much a retrieved parame- 
ter will respond to an actual change in that parameter 
in the atmosphere. For a given change in the true atmo- 
spheric state of some parameter i, if the measurements 
are perfect, we would expect to retrieve exactly that same 
change and hence, the value of An would be one. If the 
measurements contribute nothing to our knowledge of 
parameter i, that is all of our knowledge of its value is 
from the prior, then An will be zero. We can use this 
diagnostic to assess how heavily our error estimations 
are informed by the measurements. This is most impor- 
tant when using the Level-by-Level temperature profile 
retrieval. The sum of the diagonal elements of this ma- 
trix determines the total number of independent pieces 
of information that can be retrieved from the measure- 
ments. This can never be larger than the total number 
of individual data points. 

3. TEST ON SYNTHETIC MEASUREMENTS 

In this section, we apply the CHIMERA to a set of 
synthetic measurements in order to assess the robustness 
of each retrieval algorithm. 

3.1. Synthetic Observations 
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We create a generic hydrogen dominated hot jupiter 
planet and derive its emission spectrum in three different 
observing scenarios. Table IT] summarizes the basic planet 
parameters used to generate the model atmosphere and 
contrast spectrum. For simplicity we assume that the 
trace species have mixing ratios that are constant with 
altitude. Equations [I3p6] are used to generate the atmo- 
spheric temperature pronle of the planet from the val- 
ues in Table [I] . Figure [2| shows the model atmosphere 
and spectrum of the synthetic planet. The raw flux is 
divided by a PHOENIX stellar grid model that closely 
matches the chosen stellar properties. The thermal emis- 
sion contribution functions (Figure l2J top-right) suggest 
that most of the emission originatesbetween a few bars 
and a few mbars. Our synthetic data only provide be- 
lievable estimates for the temperatures and abundances 
over this region of the atmosphere. The thermal contri- 
bution functions indicate that the emission from shorter 
wavelengths comes from deeper layers in the atmosphere 
and regions of high opacity tend to push the emission to 
higher altitudes. In this example, water is the dominant 
opacity source and acts almost like a continuum absorber 
across the spectrum (Figure l2| bottom- right ). If we had 
no absorbing molecular species other than 112/He most 
emission would originate from the ~10 bar level. 

We now create simulated observations for our synthetic 
planet under three different regimes. The first regime 
is a set of broadband observations through four of the 
Spitzer Infrared Array Camera (IRAC) channels at 3.6, 
4.5, 5.7, and 8 ^m (top-left, Figure [s]). This represents 
the spectral quality that is most commonly available for 
hot jupiters today. To create the synthetic observations, 
the spectrum in Figurcjllis first integrated over the IRAC 
filter functions at each channel and then random noise is 
added to each channel determined by the error bars size. 
The size of the error bars are representative of typical 
errors of IRAC observations (e.g., Machalek et al. 2009). 

The second observational scenario is a multi- 
instrument case combining both Spitzer photometry, 
ground based photometry, and Hubble Wide Field Cam- 
era 3 spectra (WFC3) (top- right Figure l3|. This com- 
bined set of observations from various instruments is 
more representative of the current level of observations 
that can be made today, and likely for the next half- 
decade, for many planets (e.g., WASP12b, WASP4b, 
HD209458b). Again, we use the same four Spitzer IRAC 
channels and error bars as before but also include ground 
based H and K^ band photometry points. The error bars 
are taken from Crossfield et al. (2012). To create the 
synthetic WFC3 measurements (1.15-1.63 /xm), the high 
resolution spectrum is convolved with a Gaussian instru- 
mental profile with a FWHM of 0.0325 /im with error 
bars taken from Swain et al. (2012). Random noise is 
again added to each point. 

The third observational scenario illustrates the perfor- 
mance of a potential modest (by modest we mean rea- 
sonable cost) future spaceborne, FINESSE- like, telescope 
might produce (bottom. Figure J3]). These simulated ob- 
servations are created by convolving the high resolution 
spectrum with a moderate resolution Gaussian instru- 
mental profile with a FWHM of 0.0075 /xm (R~300 at 
2 /xm). The measurement error bars, and hence random 
noise, are only suggestive and are based on a FINESSE- 
like noise model (Swain 2012). This spectral resolution is 



comparable to that of Exoplanet Characterization Obser- 
vatory (Echo) below 5/i, but less than the James Webb 
Space Telescope (JWST) Near Infrared Spectrometer. 

Aside from the potential development of a ground- 
based near-IR spectroscopy program, most observations 
for the foreseeable future are likely to fall somewhere be- 
tween the first and second cases. We are also being opti- 
mistic in our "worst" case observational scenario by in- 
cluding four broadband points instead of the now typical 
two from IRAC. In the latter case it is impossible to pro- 
vide any unique constraints on the atmosphere without 
imposing many pre-conceived assumptions and priors. 

3.2. The Prior 



As mentioned in S 2.1.1 the prior is important when 
the spectral information content is limited. We use a 
prior on both the gases and the temperature profile. For 
the purposes of this synthetic study, we assume Gaussian 
priors on the parameters that control the temperature 
profile and on the gas abundances. We could have cho- 
sen flat (un-informed) priors, however, the formalism of 
optimal estimation requires that the prior be Gaussian, 
and hence we maintain this prior for all of the retrieval 
approaches. We choose extremely broad Gaussian pri- 
ors in order to mitigate the influence they have on the 
retrievals. For the temperature profile prior, we choose 
parameters that would reasonably match an atmosphere 
that is in radiative equilibrium over a wide range of con- 
ditions (e.g., variations in Km, jy-^, 7-U2, /?, and a ). Table 
[2] lists the prior parameters we use in terms of the prior 
mean, Xa and the prior covariance matrix, Sa. 

In addition to the Gaussian priors, we impose a lower 
limit on mixing ratios with a value 1 x 10^^^ and an upper 
limit requiring the sum of the mixing ratios of the four 
retrieved gases to be less than 0.15. These limits attempt 
to bound what can be reasonably expected for the com- 
positions of a hydrogen dominated atmosphere. Also, it 
would be impossible to detect a gas with an abundance 
less than 1 ppt in these simulations. We also impose 
a limit on the parameters that govern the temperature 
structure. We do not allow kjj^ to go above or below 10 
and 1 X 10"'*, respectively. The lower limit is roughly 
the order of magnitude value of the Planck mean opac- 
ity expected for an all hydrogen atmosphere. The up- 
per limit is a bit extreme but would be representative of 
an extremely opaque atmosphere. The upper and lower 
bounds on 71 and 72 are between 10 and 1 x 10~ and 
are chosen to allow for a reasonable span of tempera- 
ture profiles ranging from ones with inversions to ones 
nearly transparent to solar radiation, a can only have 
physically meaningful values between and 1. /3 can- 
not have values below 0, and we impose an artificial up- 
per limit of 2. Generally these upper and lower limits 
rarely matter as most of the posteriors lie well within 
their ranges. Figure |4] shows the resulting temperature 
distribution and gas priors (inset). The prior tempera- 
ture profile distributions are reconstructed by propagat- 
ing the Gaussian prior probability distributions (includ- 
ing the above limit s) of k j_r, 'fvi, 7u2: /^j ^^nd a in Tablel2] 
through equations |13|16[ Upon reconstructing the tem- 
perature profiles there are thousands of temperatures for 
each pressure level. With these profiles a histogram of 
temperatures at each level can be constructed. Rather 
than show the "spaghetti diagram" with thousands of in- 
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dividual profiles, we show the la (68%) and 2 a (95 %) 
confidence intervals at each pressure level. These confi- 
dence intervals are what is shown in Figure [4] 

3.3. Results from the Parameterized Temperature Profile 

We apply the three retrieval techniques to the three 
synthetic observations in Figure [3] under the radiative 
equilibrium temperature profile parameterization. This 
temperature profile has few parameters (5) because of 
the assumption of radiative equilibrium. The tempera- 
ture prior in Figure |4] is used in all three techniques for 
all three retrieval cases. In each case, we retrieve the 
parameter distributions for the following state vector: 

X = [log{fH2o)Jog{fcHi),log{fco),log{fco2), 

log{KiH),log{jy,),log{'yy^),/3,af. (18) 

where again, the fiS are the altitude independent volume- 
mixing ratios. We start by first applying the optimal 
estimation approach. In order to ensure that the retrieval 
does not get stuck in a local minimum, multiple starting 
guesses are used. These typically all converge to th e same 
temperature and gas solution. As described in ^2.1.3 



the covariance matrix and the best fit from OF are then 
used to initialize the DFMC chains. Finally, the best fit 
from OF is used to initialize the synthetic measurement 
realizations used in the BMC. figures [5][8] and Table |3] 
summarize the retrieval results and form the basis for 
the comparisons. The bounds quoted in Table [3] are for 
the 68% confidence intervals. 

We must be careful in interpreting the confidence in- 
terval values when the posteriors extend to the imposed 
upper and lower limits, especially when those limits are 
somewhat arbitrary. Parameters with posteriors that ap- 
proach the imposed lower limit will result in an overesti- 
mate of the lower bound on the confidence interval and 
an underestimate in the upper bound due to the imposed 
upper limit. In some cases if there were no Gaussian prior 
or no lower limit, the lower bound could extend to — oo! 
Of course we would interpret such a case as only having 
an upper bound. 

Figure [5 shows the spectral fits. The first row shows 
the singlebest fit from optimal estimation. The second 
and third rows show the fits from the BMC and DFMC, 
respectively. Since the BMC and the DFMC provide 
many thousands of spectra, rather than plot each one, 
we summarize the fits by showing the median spectrum 
along with the 1- and 2(7 spread at each wavelength. Il- 
lustrating the fits in this manner is more representative 
of the posterior than plotting spectra of different chi- 
squared levels. In other words, if a random set of pa- 
rameters is drawn from the posterior, there would be a 
95% chance that the flux at any one wavelength of the 
spectrum resulting from that parameter draw would fall 
with in the 2a spread etc. There is little if any spread in 
the fits as the measurement quality improves. In the fol- 
lowing subsections we summarize posteriors for the gas 
compositions, temperatures, and C to O ratios for each 
observational scenario. 

3.3.1. Gas Abundance Retrievals 

Figures |6] and [7] summarize the gas mixing ratio 
retrieval results. Figure [6] shows histograms of the 
marginalized posterior for each of the four gases as a 



result of each retrieval approach and observational sce- 
nario. We take the DFMC posterior (blue) to be repre- 
sentative of the true posterior. The optimal estimation 
posteriors (solid red curve) and the prior (dot-dashed red 
curve) are smooth because they are constructed analyti- 
cally from the diagonal elements of S. Figure[7|is a "stair- 
step" plot that shows the correlations amongst the four 
gases comparing optimal estimation to the DFMC. For 
brevity, we do not show the Bootstrap Monte Carlo corre- 
lations. The solid blue filled regions are the 1 (dark) and 
2 (light) a confidence intervals derived from the DFMC, 
and the red curves are the la (inner) and 2cr(outer) con- 
fidence intervals derived from OF. 

The first column of Figure [6] shows the marginalized 
gas posteriors for the broadband observational scenario 
and the top set of panels of Figure [7] show the correla- 
tions amongst the gases. In this observational scenario, 
the three retrieval techniques produce quite different pos- 
teriors. H2O has a fairly narrow posterior (relative to the 
prior) near the true state, suggesting that it is reasonably 
well constrained, at the la level, by even this low infor- 
mation content spectra. This is unsurprising as water 
is prevalent across all of the channels in this spectrum 
(see Jacobian in Figure [2]) . At the 2a level, however, 
the OF retrievals provide less of a constraint (Figure [7 1 . 
The CH4 and CO2 posteriors abruptly fall off towards 
their upper end suggesting an upper bound constraint 
on these gases. The low end of their posteriors begin 
to track the prior down to the imposed lower limit indi- 
cating that from this observational scenario, there really 
is no observable lower limit to the abundances of these 
species. CH4 has a better defined upper edge than CO2 
because both the 3.6 ^m and 8 /im channels overlap with 
the strongest methane absorption bands. CO is virtually 
unconstrained by the synthetic broadband measurements 
as it closely matches the prior across the full range of val- 
ues. The difficulty in constraining CO and CO2 is due to 
the inability of the 4.5 fim broadband photometric mea- 
surement to decouple the CO and CO2 strong bands. 
There is a slight hint of an inverse correlation, as ex- 
pected, from the optimal estimation results in Figure [Tj 
however this correlation is not apparent from the DFMC 
results. 

Our implementation of optimal estimation struggles to 
appropriately capture the errors in this observational sce- 
nario. This is because it approximates the posteriors 
with broad Gaussians which simply do not capture the 
appropriate structure. It does, however, do a fairly good 
job of determining the true state. The Gaussian approx- 
imation cannot appropriately handle upper bounds on 
CII4 and CO2, causing an overestimate of the la up- 
per bound. We note, that at least it is overestimating 
the errors rather than underestimating them. OF does 
a fine job at approximating the posterior for CO, which 
happens to be similar to the prior, again suggesting no 
constraint. This is reaffirmed by looking at the aver- 
aging kernel element. A, in Table [3] which shows that 
most ( 70%) of the i nforma tion in constraining CO comes 
from the prior (see ^2.2.1 1. The Bootstrap Monte Carlo 
(green curve) completely fails to appropriately capture 
the posterior in this particular observational scenario. 
This is because the different realizations from the BMC 
approach produce parameter distributions that are con- 
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fined to a small area of phase space localized near the best 
fit solution from optimal estimation and cannot therefore 
sample the entirety of the posterior. 

The second column of Figure [6] shows the gas posteri- 
ors for the multi-instrument observational scenario. The 
information gain from this synthetic observational sce- 
nario is only marginal relative to the broadband case. 
Water has the largest improvement in uncertainty due to 
the leverage provided by the WFC3 spectra which cov- 
ers the 1.15 fi and 1.4 /im water bands. Upon inspecting 
Figure J7] we find that the WFC3 data combined with the 
groundDased photometric points trims the 2cr tail but 
does little to improve the la uncertainties. Sadly, there 
is virtually no reduction in the uncertainties on CH4, 
CO, and CO2. In fact the marginalized posteriors (Fig- 
ure pi) produced by DEMC in this observational scenario 
look nearly identical to the previous case. These results 
are not surprising given that these species do not have 
strong features at the near-IR wavelengths of the WFC3 
spectroscopy and the H and K^ photometry. Optimal es- 
timation provides an accurate error estimation for water 
but appears to provide an overly optimistic estimation 
of the uncertainties on CH4. OF, as in the previous sce- 
nario, captures the essence of large uncertainties on CO 
and CO2 with broad Gaussians. Bootstrap Monte Carlo 
underestimates the uncertainties in all species with the 
exception of water. 

Finally, results for a hypothetical future spaceborne 
telescope are shown in the last column of Figure |6] and 
bottom set of panels of Figure [7] The reduction in un- 
certainties are staggering when compared to the previous 
observational scenarios. All of the gases are constrained 
to within better than an order of magnitude. This is 
an orders-of-magnitude improvement over the previous 
cases. The high signal-to-noise and high spectral resolu- 
tion combine to provide excellent coverage of each of the 
four gases. With this high quality spectrum, all three 
retrieval approaches give the same results. The Gaus- 
sian assumption used in optimal estimation is perfectly 
appropriate in this case. The differences in the la uncer- 
tainties derived from OE are less than ~ 10% than the 
uncertainties derived from DEMC. The two-dimensional 
confidence intervals in FigurelTJalso agree quite well. The 
prior also plays very little role in the retrieval as shown 
with the near unit averaging kernel elements in Table [Sj 

3.3.2. Temperature Profile Retrievals 

Figure [8] shows the temperature profile retrieval results 
{la (dark red) and 2a (light red)) reconstructed from the 
posteriors of log{Kiii),log{'~^^^),log('^y^), f , and a, as in 
Figure |4] The blue curve is the statistical median of 
these profiles. This median profile not representative of 
any one given temperature profile, and in fact this me- 
dian profile may not even provide a good fit to the obser- 
vations or adhere to the parameterization, but is shown 
simply as a statistical summary of all possible temper- 
ature profiles. The black curve is the true temperature 
profile. 

The temperature profile posteriors for the broadband 
scenario (first column Figure Isl) from OE and DEMC 
have similar widths and both capture the entire true 
temperature profile within the 2a interval. There is also 
a non-negligible (^30-50%) reduction in the tempera- 
ture precision compared with the prior, over the atmo- 



spheric region probed by these observations. Outside 
of the range spanned by the thermal emission contribu- 
tion functions (Figure 2]), the temperature uncertainty 
grows and begins to relax back towards the prior where 
there is no observational constraint. Again, the BMC ap- 
proach completely underestimates the error when com- 
pared with the other two approaches because of its inabil- 
ity to fully characterize the posterior outside of a small 
region of phase space localized around the OE original 
best fit. 

Moving onto the multi-instrument observations (mid- 
dle column Figure Is]) we find a 22% reduction in the tem- 
perature uncertainty between 1 and 0.01 bars. The large 
number of spectral channels from the WFC3 data that 
have weighting functions over this region are the primary 
contributors to this increased precision. OE and BMC 
underestimate the temperature uncertainties relative to 
DEMC at 100 mbars, but the OE and DEMC have rea- 
sonable agreement over the entire profile. 

The future spaceborne telescope observations improve 
the temperature uncertainties by a remarkable factor of 
^ 4.5 over the previous case. Optimal estimation slightly 
overestimates the temperatures outside the atmospheric 
levels probed by the observations. This overestimate is 
due to the overestimation of the 7„i and 7^2 posterior 
widths. These two parameters control the relative dif- 
ference between the upper atmosphere and lower atmo- 
sphere temperatures. Hence, extreme values of 7j,i and 
7t,2 will effect these regions more than the middle atmo- 
sphere. This why the OE and the DEMC temperatures 
agree in the middle atmosphere but not outside of it. 

3.3.3. C to O Ratios 

Determination of the C to O ratios of explanatory at- 
mospheres is critical to the understanding of their at- 
mospheric chemistry (Lodders & Fegley 2002; Moses et 
al. 2012) and formation environments (Oberg, Murray- 
Clay, & Bergin 2011). Given the abundance posteriors 
derived with CHIMERA, we can compute C/0 posterior 
distributions (Figure[9]) . The C/0 is calculated with the 
following formula, 



C/O 



EC _ CHi + CO + CO2 

Yn ^ H2O + C0 + 2CO2 ' 



(19) 



There are a few simple things to note about this equa- 
tion. When CO is the dominant species, C/0 is 1. If 
CO2 is the dominant species, C/0 will be 1/2. When 
methane dominates the C/0 will be large and when H2O 
dominates C/0 will be small. The solar C/O is 0.55. A 
number of exoplanets have reported C/O's near 1. We 
can construct the C/0 probability distributions by prop- 
agating t he p robability distributions of each gas through 
equation |19[ similar to the method used to construct 
the temperature profile posteriors. Before inspecting the 
posteriors derived from CHIMERA, we find it illustrative 
to investigate the prior. Upon propagating the Gaussian 



priors (with the limits) of the gases through equation 19 
we obtain the C/0 prior in Figure l9] Wc find that this 
prior has two peaks, one is at a C/O of 1 and the other 
is at a C/O of 0.5. The locations of these peaks are in- 
sensitive to whether or not the gas abundance priors are 
uniform or broad Gaussians. These peaks are also insen- 
sitive to the lower and upper bounds placed on either a 
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uniform or Gaussian gas prior. These double peaks are 
due to an elegant mathematical misfortune. For illustra- 
tive purposes, let us assume we draw the set of four gases 
from a uniform log prior. We would expect then, that 
one gas will have a larger abundance than the other three 
25% of the time. This means that in roughly 25% of the 
draws CO will dominate, which would cause the ratio in 



equation 19 to be one, 25% of the time. A similar argu- 
ment goesfor CO2 resulting in a C/0 of 0.5 roughly 25% 
of the time. So we see, that if we did not observe a par- 
ticular planet, and assumed uniformed priors on each of 
the gases, we would naturally conclude that the planet 
has equal chances of having a C/0 of one or one-half. 
However, our priors are not uniform, but rather broad 
Gaussians, but they are broad enough that this behavior 
still occurs, with a slight preference for a C/0 of one. 

If we observe this planet with the four broadband 
points we obtain the posteriors in the upper-left panel 
of Figure [9] OE and DEMC produce similar C/0 poste- 
riors both of which maintain the two peaked features at 
1 and 0.5 but with overall less power in the peaks. These 
features persist simply because the gas posteriors from 
DEMC and OE do not deviate too strongly from the 
prior. There is more probability in the lower C/0 tail 
than in the prior because of the higher values of H2O 
preferred by the measurements over the prior. This is 
good, since the C/0 for our fictitious planet is much 
less than one or one-half. BMC anomalously captures 
the true C/0 at the peak of its posterior. Again, BMC 
greatly underestimates the posterior widths because it 
only searches a localized region about the OE best fit. 
Since the OE best fit gas abundances are very near truth 
(Figure IgI), the BMC posteriors, which are highly local- 
ized about the OE best fit parameters, will overempha- 
size the C/0 derived by that best fit. 

The story is the same for the multi-instrument obser- 
vational scenario. Unfortunately, at least in this exam- 
ple, it appears that the WFC3 and ground based data 
provide very little additional constraints in reducing the 
C/0 uncertainty, with the double peaked feature from 
the prior persisting in the DEMC results. 

Improving the observational quality further with a fu- 
ture spaceborne telescope essentially obtains the correct 
value to high precision. The peak of the posterior is 
far enough away from the double peaked prior that the 
results appear to be less infiuenced by the prior than 
the previous cases. All three retrieval approaches give a 
nearly identical posterior. 

From this exercise we have learned that it is difficult to 
constrain the C to O ratio of an exoplanet atmosphere. 
Simple, uninformed, or nearly uninformed priors on the 
gas abundances produce a double peaked C/0 prior at 
near solar value and one. Even in the best cases, current 
observations are likely to provide only an upper limit on 
this quantity. This result suggests that previously pub- 
lished claims to detect enhanced C to O ratios with pho- 
tometry alone may be influenced by these subtle biases, 
and should be viewed with strong skepticism. We will 
discuss this issue in the context of specific exoplanets in 
a companion paper. 

3.4. Results from the Level-by- Level Temperature Profile 

The Level-by-Level temperature profile approach at- 
tempts to determine the temperature for each model 



layer. In contrast to the parameterized temperature pro- 
file which had only 5 parameters, the Level-by-Level ap- 
proach requires as many parameters as model layers, for a 
total of 90 parameters. This larger number of parameters 
is far greater than the number of meaningful constraints 
provided by most observations. However, this approach 
makes no assumptions about the physical structure of the 
temperature profile (e.g., radiative, radiative-convective, 
advection etc.). While there is no potentially biasing 
parameterization, the retrievals can result in unphysical 
temperature profiles. Obviously, the temperature at each 
of the 90 levels cannot be perfectly retrieved, but rather 
the retrievals have to depend on the prior when spectral 
information on the temperature is sparse. 

For the Level-by-Level prior we assume an a prior vari- 
ance of 400 K and covariance amongst each level with all 
other levels given by equation [T2j The 400 K variance 
is used to produce a similar temperature profile prior as 
in Figure |4] This correlation helps reduce the effective 
number oiTevels that have to be independently retrieved. 
Admittedly, the degree of correlation is somewhat of an 
external arbitrary parameter, but it is chosen to avoid 
over-fitting (i.e., fitting to the noise) without hindering 
the Level-by-Level flexibility. It can be thought of as a 
smoothing, or more specifically, a regularization. We can 
also use the averaging kernel profile to asses where the 
temperature is constrained by the measurements versus 
the prior. The gas priors are the same as before. We 
choose only to compare the results from optimal estima- 
tion and the Bootstrap Monte Carlo. We do not attempt 
the DEMC approach on such a large (^^100) number of 
parameters, as MCMC algorithms (to the best of our 
knowledge) are not well suited for large numbers of pa- 
rameters because of the large number of steps required 
to fully map the n-dimensional probability distribution 
when n is large. 

Figure [T0| shows the spectral fits as a result of OE and 
BMC using the Level-by-Level temperature approach, 
similar to Figure [5] Figure 11 shows the marginalized 
gas posteriors. Wefind the gas posteriors and the agree- 
ment amongst the retrieval techniques are very similar 
to those derived in Figure |6] using the parameterized 
temperature profile. This is somewhat surprising given 
the extremely different temperature profile retrieval ap- 
proaches. This suggests that the gas abundances can be 
properly and consistently retrieved regardless of the tem- 
perature profile assumptions. We could, however, imag- 
ine a contrived example in which the true temperature 
profile is so wildly different from what can be reasonably 
approximated with the parameterization, that the two 
approach es w ould yield differing gas posteriors. 

Figure 12 shows temperature profile posteriors under 
the Level-by-Level temperature profile assumption. For 
the broadband scenario, the uncertainties more or less do 
not improve much beyond the prior. The greatest gain in 
improvement is over the region spanned by the averaging 
kernel (green curve). The uncertainty reduces from the 
prior uncertainty of ±400 K to ±265 K at 100 mbars. 
The BMC approach using the Level-by-Level tempera- 
ture profile produces a much smaller error than the OE 
approach, and at some levels the 2ct uncertainties do 
not even capture the true state. Again, this because 
the BMC is only able to characterize a highly localized 
region around the OE best fit. 
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The reduction in temperature uncertainty due to the 
addition of the WFC3 and ground-based photometry 
data is more apparent with the Level-by-Level approach 
than with the parameterized approach. The uncertain- 
ties in temperature at 100 mbars are reduced to ±177 K, 
though smaller uncertainties are achieved at deeper lev- 
els due to the addition of the WFC3 data which probe 
deeper atmospheric levels. This is why the averaging 
kernel profile peaks at a deeper level. The BMC results 
show a larger uncertainty than they do in the broadband 
observational scenario but still greatly underestimate the 
profile spread relative to optimal estimation. 

The future spaceborne telescope observations reduce 
the temperature uncertainties to ±70K, a factor of nearly 
four better than what can be done with the broadband 
observations. Outside of the region spanned by the aver- 
aging kernel uncertainties relax back to the prior widths. 
As before, in both cases the BMC approach underesti- 
mates the temperature uncertainties relative to the OE 
derived uncertainties. The uncertainties in temperature 
derived using the Level-by-Level temperatures are a fac- 
tor of two larger than those derived with the parameter- 
ized temperature profile. This is because the retrievals 
with the parameterization only allow temperature pro- 
files that conform to radiative equilibrium whereas the 
Level-by-Level retrievals can allow for a wider range of 
possibilities that do not necessarily have to conform to 
this constraint. 

Another way to determine the robustness of the Level- 
by-Level retrieval is to explore the role of the prior tem- 
perature profile (see e.g., Lee et al. 2012). For this, we 
investigate the effect of three different temperature pri- 
ors (different prior profiles, Xa but same widths, Sa) on 
the retrieved profiles and check to see if they are consis- 
tent with the estimated errors (Figure 13). The shaded 
grey region in Figure 
from Figure [12] using t 



13] is the la retrieval uncertainty 
the nominal prior. Two of the other 
priors are the nominal profile with a ±500 K offset, and 
the third is an isothermal profile set to the equilibrium 
temperature of the planet. In all three cases we find 
that the retrieved profiles fall with in the la bounds 
of the nominal retrieved profile. This suggests that al- 
though different temperature profile priors are used, they 
generally produce retrieved profiles that are statistically 
consistent with each other. As the spectral quality im- 
proves, the different priors produce identically the same 
retrieved profiles over the atmospheric regions spanned 
by the thermal emission weighting functions. Outside of 
this region, the profiles diverge and relax towards their 
respective priors with no consequence on the spectra. 
This is yet another demonstration that the spectra are 
only sensitive to a small region of the atmosphere be- 
tween a few bars up to a few mbars. While some of these 
Level-by-Level profiles may not be physical, especially in 
the broadband observational scenario, they are a more 
direct reflection of the information provided by the mea- 
surements in the absence of a parameterized model. 

4. DISCUSSION & CONCLUSIONS 

We have developed a new statistically robust suite 
of exoplanet atmospheric retrieval algorithms known as 
CHIMERA. T his su ite consists of the Optimal Esti- 
mation (OE, j]2.1.1| ), Bootstrap Monte Carlo (BMC, 



i 2.1.2), and Ditterential Evolution Markov Chain Monte 



Carlo (DEMC . j]2.1.3 l approaches and a validated for- 
ward model (j ]2.^ We have tested each of these ap- 
proaches on the dayside thermal emission spectra for a 
synthetic planet under a variety of observational scenar- 
ios ranging from current observations to potential future 
observations (@- In general, we find that the three re- 
trieval approaches produce similar posteriors when the 
measurement quality is good, typically when there are 
more observed spectral channels than retrievable param- 
eters (Figures p\ and JtI ) . The Gaussian approximation 
made by Optimal Esflmation breaks down for low qual- 
ity measurements, but is perfectly valid for high qual- 
ity measurements likely to come from future spaceborne 
observations. It also appropriately captures the correla- 
tions amongst the various parameters. This approach is 
much less of a computational burden than Monte Carlo 
approaches and will prove useful for quick reductions of 
large, high-quality data sets. The Optimal Estimation 
formalis m also allows for the calculation of the averaging 
kernel {{ 2.2.1 ), which is a useful diagnostic to determine 
how much ot the posterior is infiuenced by the prior ver- 
sus the measurements. The Bootstrap Monte Carlo ap- 
proach generally fails to capture the essence of the poste- 
riors. This is because the regeneration of synthetic mea- 
surement realizations based on the Optimal Estimation 
best fit only sample a localized region of phase space near 
the best fit solution. This is especially problematic in the 
cases where there are fewer spectral data points than pa- 
rameters. In this scenario, even with small measurement 
error, there will still be many possible best-fit solutions, 
thus creating enormous degeneracies among the parame- 
ters. Since the BMC is initialized with only one possible 
best-fit set of parameters out of many, the derived pa- 
rameter uncertainties will only be representative of the 
uncertainties about that localized best fit. We therefore 
strongly advise against the Bootstrap Monte Carlo ap- 
proach when the number of parameters is larger than 
the number of spectral data points. In the high signal- 
to-noise and high spectral resolution regime, both the 
BMC and OE methods provide reasonable parameter un- 
certainties. We have also introduced the application of 
Differential Evolution Markov Chain Monte Carlo to the 
spectral retrieval problem and found that convergence 
can be obtained efficiently by using an appropriate pro- 
posal distribution based on the chain history. 

We find that for broadband observations typical of 
most hot jupiters, very limited constraints on the gas 
abundances can be obtained. The Spitzer photometry 
does a particularly poor job constraining the relative 
abundance o f CO, with most posteriors simply reflecting 
their priors (^3.3.1 ). The addition of WFC3 observations 
provide little additional constraint on the gas abundances 
derived from dayside thermal emission spectra, with the 
exception of a slight improvement on the water abun- 
dances. This is primarily due to the limited spectral cov- 
erage provided by the red grism on WFC3, which spans 
the wavelengths from 1.2-1.6 /xm. As the measurement 
quality improves, the parameter uncertainties decrease 
and become more Gaussian. Moderate cost future space- 
borne instruments have the capability of obtaining better 
than order-of-magnitude constraints on gas compositions 
with their posteriors generally being independent of the 
prior. This is typically many orders of magnitude better 
than current observational capabilities (Table |3| . The 
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derived gas posteriors are also independent of whether 
or not a parameterized or Level- by-L evel temperature 
profile is used (Figure [6 vs. Figure 111. 

Constraining the C to O ratio of exoplanet atmo- 
spheres is very difficult due to the broad nature of some 
of the gas posteriors, especially CO. In the absence of 
valid observational constraints the posteriors for these 
molecules simply reflect the priors, which produce a 
double-peaked disti ributi on with maxima at C to O ra- 
tios of 0.5 and 1 (j ]3.3.3 Figure [o]). Only high quality 
observations from the future spaceborne telescope scene- 
rio are independent of the double-peaked prior. As a 
result, caution must be taken when interpreting C to O 
ratios from broad gas posteriors. 

Reasonable temperature constraints could be obtained 
in all observational scenarios and temperature retrieval 
approaches, though the Bootstrap Monte C arlo ap proach 
again fails to fully capture the posterior ([3.3.2 Figure 
^. The temperature profiles and corresponding uncer- 
tainties can only be trusted for the region over which the 
thermal emission contribution functions peak, typically 
between a few bars and a few mbars (Figure l2| . Out- 
side of this window, the temperature profiles are strongly 
affected by their priors. The Level-by-Level tempera- 
ture profile approach overestimates the temperature un- 
certainties compared with the p aram eterization due to 
the allowance of more profiles (j |3.4[ ). These Level-by- 
Level profiles can be unphysical but are more reflective 
of the measurements without imposing preconceived no- 
tions of how the physical structure of the atmosphere 
should behave. While this approach produces statisti- 



cally consistent profiles in low quality observational sce- 
narios, we would still recommend using a parameteriza- 
tion for said cases. However, for high quality spectra 
the Level-by-Level approach is recommended given its 
slightly more pessimistic temperature uncertainties and 
its non-dependence on a particular parameterization. 

In a follow up investigation we will use CHIMERA to 
perform a uniform analysis of an ensemble of secondary 
eclipse spectra. Such a study will allow us to determine 
the biases introduced by the choice of fitting method for 
individual planets and to derive a uniform set of relative 
abundances and temperatures for these planets that can 
be reliably inter-compared and trends identified. This 
kind of uniform analysis has the potential to provide 
invaluable insights into exoplanetary atmospheric pro- 
cesses and formation environments. 
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TABLE 1 
Parameters used to generate the fictitious model atmosphere and spectrum. Rp is the planet radius in units of Jupiter 

RADII, Rstar IS THE STELLAR RADIUS IN UNITS OF SOLAR RADII, Tstar IS THE STELLAR EFFECTIVE TEMPERATURE, a IS THE SEMIMAJOR AXIS, 

Tint IS THE INTERNAL HEAT FLUX TEMPERATURE OF THE PLANET, g IS THE PLANETARY SURFACE GRAVITY. 7i,i, 7„2, ^^IR, a, AND /3 ARE 

THE PARAMETERS THAT CONTROL THE SHAPE OF THE RADIATIVE EQUILIBRIUM TEMPERATURE PROFILE. THE fi'S ARE THE 

CONSTANT- WITH-ALTITUDE VOLUME MIXING RATIOS FOR EACH SPECIES IN PARTS-PER-MILLION (pPM). 



Parameter 


Value 


Rp{Rj) 


1.138 


J^star y-ttsun ) 


0.756 


Tstar (K) 


5040 


a{AU) 


0.031 


Ti„t{K) 


100 


log{g){cms 2) 


3.341 


Ivl 


1.58 X 10-1 


1v2 


1.58 X 10-1 


Kin{am?s-^) 


3 X 10-2 


a 


0.5 





1.0 


fH2 (ppm) 


8.5 X 105 


fHe (ppm) 


1.5 X 10^ 


fH20 (ppm) 


370 


fcHi (ppm) 


1 


fco (ppm) 


31.6 


fc02 (ppm) 


0.2 



TABLE 2 
Gaussian prior parameter values and widths. The true state is the same as in Table [l] but in logarithmic units for some of 

THE parameters. ThE MIXING RATIOS OF EACH GAS, fk, ARE IN PPM. THE INFRARED OPACITY, Kjji, HAS UNITS OF cm?g~^. 7„i, 7„2, O, 
AND /3 ARE ALL UNIT-LESS. NOTE THAT WE RETRIEVE THE LOG OF ALL VALUES EXCEPT a AND /3. 



Parameter 

'09(7^1) 

log("iv2) 

loginin) 

a 

/3 

log{fH2o) 

logifcHi) 

logifco) 

I09UC02) 



True 
-0.8 
-0.8 

-1.52 

0.5 

1 

2.568 
0.0 

2.663 

-0.70 



Prior State (xa^i) 
^09 
-0.7 
-2.0 
0.5 

1 

2 

1 

2 

1 



Prior Width U{Sa,u)) 

1 

1 

0.5 
0.05 
0.25 

6 

6 

6 

6 
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TABLE 3 

Numerical summary of the retrieval results for several parameters as derived from each retrieval technique and 

observational scenario. for each parameter and each observational scenario we show the true value, the la (68% 

confidence interval) marginalized prior uncertainties, and the 1(t marginalized uncertainties derived from optimal 

Estimation (OE), Bootstrap Monte Carlo (BMC), and Differential Evolution Markov Chain Monte Carlo (DEMC) as 

WELL AS THE AVERAGING KERNAL ELEMENT FOR THAT PARAMTER (A) . THE GAS ABUNDANCES, fi, ARE GIVEN IN TERMS OF VOLUME MIXING 

RATIO. We also show a REPRESENTATIVE TEMPERATURE (100 MBARS TEMPERATURE) AND THE C-TO-O RATIO. THIS TABLE IS LAYED OUT 

SO THAT FOR A GIVEN PARAMETER EASY COMPARISONS IN EITHER THE OBSERVATIONAL SCENARIO (lEFT-RIGHT) OR THE RETRIEVAL 

TECHNIQUES (tOP-BOTTOm) CAN BE MADE. 



Parainter 



Broadband 



Multi-Instrument 



3.70 X 10""-* 
10-10 -3.92 X 

06 _ 7.82 X 
05 _ 2.68 X 

07 _ 1.24 X 
0.872 



3.70 X 10-04 
4.94 X 10-10 _ 3.92 x IQ-o^ 
3.68 X 10-06 - 7.03 x lO-""* 
1.18 X 10-05 -2.45 X 10-04 
7.06 X 10-06 - 1.40 X 10-0'^ 
0.983 



Future Telescope 



3.70 X 10-04 
4.94 X 10-10 _ 3.92 x IQ-oa 
2.31 X 10-04-4.56 X 10-04 
2.44 X 10-04 - 4.03 X 10-04 
2.83 X 10-04 - 4.98 x 10-04 
0.999 



/. 



H20 



True; 

Prior: 

OE: 

BMC: 

DEMC: 

A: 



4.94 X 
1.25 X 
5.74 X 
1.84 X 



10- 
10- 
10- 



10-03 
10-03 
10-04 
10-03 



fcH4 



fco 



True: 

Prior: 

OE: 

BMC: 

DEMC: 

_A^ 

True: 

Prior: 

OE: 

BMC: 

DEMC: 

A: 



2.63 X 
8.89 X 
6.40 X 
2.21 X 



4.94 X 
4.97 X 
2.05 X 
2.00 X 



1.00 X 10-06 
10-10 -2.55 X 
10-11 -5.73 X 
10-0** - 1.64 X 
10-11 -9.27X 

0.259 
3.16 X 10-05 
10-10 -3.92 X 
10-10 -2.43 X 
10-06 - 6.10 X 
10-10 - 5.32 X 

0.316 



10-03 
10-03 
10-05 
10-04 



1.00 X 10-06 
2.63 X 10-10 - 2.55 x 10-03 
4.54 X 10-09 - 3.92 x 10-05 
9.07 X 10-08 - 3.89 x 10-06 
4.40 X 10-11 - 1.96 X 10-06 
0.979 

3.16 X 10-05 
4.94 X 10-10 - 3.92 x 10-03 
7.40 X 10-10 - 4.13 X 10-03 
3.00 X 10-07 - 1.81 X 10-05 
8.89 X 10-11 -7.16 X 10-05 
0.176 



1.00 X 10-06 
2.63 X 10-10 - 2.55 x 10-03 
7.50 X 10-07 - 1.54 X 10-06 
7.62 X 10-07 - 1.46 X 10-06 
7.39 X 10-07 - 1.52 X 10-06 

0.999 

3.16 X 10-05 
4.94 X 10-10 - 3.92 x 10-03 
5.87 X 10-06 - 3.13 X 10-05 
4.37 X 10-06 - 2.76 x 10-05 



3.52 X 10- 



2.66 X 10- 



0.996 



2.00 X 10-07 
2.63 X 10-10 - 2.55 x 10-03 
1.94 X 10-07 - 4.82 X 10-07 
2.03 X 10-07 - 4.29 X 10-07 
2.30 X 10-07- 5.06 X 10-07 
0.999 



fc02 



True: 

Prior: 

OE: 

BMC: 

DEMC: 

A: 



2.63 X 
7.73 X 
9.44 X 
2.21 X 



2.00 X 10-07 
10-10 -2.55 X 
10-11 - 1.78 X 
10-00 - 3.61 X 
10-11 -9.27X 
0.508 



10-03 
10-04 
10-07 

10-07 



2.00 X 10-07 

2.63 X 10-10 - 2.55 x 10-03 
5.07 X 10-00 - 4.25 X 10-03 

1.64 X 10-00 - 1.10 X 10-07 
2.33 X 10-11 -7.37X 10-07 

0.689 



Tioon^blK] True: 
Prior: 
OE: 
BMC: 
DEMC: 



1313 
876 - 1503 
932 - 1358 
1150 - 1249 
1050 - 1355 



1313 
876 - 1503 
1075 - 1274 
1117-1284 
1136 - 1374 



1313 
876 - 1503 
1267 - 1340 
1278 - 1327 
1294 - 1348 



8.00 X 10-^ 
3.82 X 10-02 - 8.00 
6.26 X 10-03 - 0.938 
1.29 X 10-02 - 0.427 
1.93 X 10-04 - 0.721 



8.00 X 10-2 

3.82 X 10-02 - 8.00 

2.25 X 10-2 - 8.93 x 10-2 

1.81 X 10-02 - 8.09 X 10-2 

1.27 X 10-02 -6.90 X 10-2 



C/O 



True: 

Prior: 

OE: 

BMC: 

DEMC: 



8.00 X 10- 
82 X 10-02 - 8. 
07 X 10-03 - 1 
60 X 10-02 - 
53 X 10-04 - 



00 
00 
32 
97 
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Fig. 1. — Comparison of the thermal emission spectrum from our forward model (black) with the NEMESIS forward model (red). The 
temperature-pressure profile is shown in the inset. For this comparison we assume uniform mixing ratios of lO"** for CH4 , CO2, CO, and 
H2O. H2 is set to 0.85 and He is set to 0.15. This planet is assumed to be hydrogen dominated (mean molecular weight of 2.3 amu) with 
a radius of IRj, a gravity of 22 ms~^ orbiting a 5700 K pure blackbody star with a radius of IRsun. 
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Fig. 2. — Synthetic planet atmosphere and spectrum. Top-Left; Model temperature-pressure profile. The solid curve is the temperature 
profile and the dashed curve is the averaged thermal emi ssion co ntribution function, or where the emission in the atmosphere is coming from. 
The temperature profile is constructed using equations |13|16] and the parameters in Table ^ Top- Right; Thermal emission contribution 
function. This plot shows where the emission is coming from as a function of wavelength, smoothed to a resolution of 0.05 microns. Red 
corresponds to the peak of the thermal emission weighting functions, where the optical depth is unity, and blue represents zero emission. 
Most emission emanates between a few bars and 0.01 bars with deeper layers probed by shorter wavelengths. Bottom-Left; Resulting 
spectrum smoothed to a resolution of 0.05 microns. Blackbodies for the hottest, coolest, and average temperatures are shown. The dotted 
curves at the bottom are the filter profiles for typical photometric observations. Bottom-Right; Gas Jacobian generate from equation 
[9] This plot shows the sensitivity of the flux contrast as a function of wavelength to the various absorbers (the units are arbitrary but 
consistent). 
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Fig. 3. — The spectrum of the synthetic hot jupiter observed in three different scenarios. These "observations" are created by convolving 
the high resolution spectrum in Figure [T] with the appropriate instrumental profiles. Random noise is then added to each data point. 
Top-Left: Synthetic observations as viewed through the Spitzer broadband 3.6, 4.5, 5.7, and 8 fim channels. Top-Right: Multi instrument 
observations that include WFC3 (1.15-1.63 ^m), ground based H and Ks, and Spitzer Broadband (3.6, 4.5, 5.7, and 8 fim). Bottom: 
Hypothetical future spaceborne observations. The dotted curves on the bottom of each plot are the photometric transmission functions. 
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Fig. 4. — Temperature and gas priors (inset). Dark red represents the la spread in the allowed temperature profiles as a result of the 
prior parameter distributions in Table [2] Light red is the 2cr spread allowed in the temperatureprofiles. The blue curve is the median 
temperature profile and the black curve is the temperature profile constructed from Xa in Table p] . The gas priors are broad Gaussians. 
H2O and CO have the same prior mean, CH4 and CO2 have the same prior mean. Note that the prior is gaussian in log of the mixing 
ratios. The y-axis on the inset is normalized probability. 
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Fig. 5. — Fits to the three different sets of data (columns) from each of the three different retrieval techniques (rows). The first scenario 
consists of the four IRAC photometry channels. The second scenario consists of the four IRAC photometry channels, ground based H and 
Ks band photometry, and HST WFC3 spectroscopy. The third scenario is representative of a FINESSE-like future, spaceborne telescope. 
For the Optimal Estimation, the single best fit is shown in blue. The Bootstrap Monte Carlo and the Differential Evolution Markov Chain 
Monte Carlo approaches generate many thousands of spectra. The median of these spectra is shown in blue and the 1- and 2a spread in the 
spectra are shown in dark and light-red, respectively. The dotted curves at the bottom of each panel are the broadband filter transmission 
functions. The insets are a zoom in of a spectral region between land 2 fim to better show the spread in the spectra. Note that there is 
virtually no spread in the spectra for the future telescope case. 
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Fig. 6. — Marginalized posterior probability distributions for each of the retrieved gases (rows) and observational scenario (columns). 
In each panel the probability distribution for each retrieval technique are shown in different colors. The Gaussian probability distributions 
from Optimal Estimation are in red, Differential Evolution Markov Chain Monte Carlo in blue, and Bootstrap Monte Carlo in green. The 
priors for each gas are the dot dashed red curve. The true answer is the vertical black line. 
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Fig. 7. — Gas correlations for each of the observing scenarios. The red curves in each are the analytic confidence intervals from the 
Optimal Estimation posterior covariance matrix (S). The inner ellipses are the la (68%) and the outer ellipses are the 2a (95%) confidence 
interval. The 1- and 2a confidence intervals derived from the Differential Evolution Markov Chain Monte Carlo are shown in dark and light 
blue, respectively. Note that the scales for the confidence intervals derived from the broadband observations (top) and the multi-instrument 
observations (middle) are the same. The scale on the future telescope (bottom) confidence intervals is much smaller. 
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Fig. 8. — Temperature profile posteriors for cacfi observational scenario (colu mns) a nd each retrieval technique (rows). The solid black 
curve in each panel is the true temperature profile constructed with equations |13|16| and the parameters in Table IT] The dashed black 
curve is constructed from the temperature parameters, Xa, just as in Figure 111 The blue curve is the median temperature profile. The dark 
and light red regions are the 1- and 2a (68% and 95%) uncertainties, respectively. 
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Fig. 9. — C to O ratio posteriors. The dot-dashed green curve is the prior, the sohd red curve is from OE, blue is from DEMC, and green 
is from BMC. The vertical dashed line is the true C/O. The top left panel are the C/O's derived from the broadband observational scenario, 
top right are the C/O's derived from the multi-instrument scenario, and the bottom are the C/O's derived from the future spaceborne 
telescope scenario. Though it appears that the BMC characterizes the C/O errors well, it is for the wrong reasons. See i]3.3.3| 
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Fig. 10. — Fits to the three different sets of data (columns) from two of the retrieval techniques (rows) using the Level-by-Level 
temperature approach. For the optimal estimation approach, the single best fit is shown in blue. The Bootstrap Monte Carlo approach 
generates many thousands of spectra. The median of these spectra is shown in blue and the 1- and 2a spread in the spectra are shown in 
dark and light-red, respectively. The dotted curves at the bottom of each panel are the broadband filter transmission functions. 
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Fig. 11. — Marginalized posterior probability distributions for each of the retrieved gases (rows) and observational scenario (columns) 
using the Level-by-Level temperature profile approach. In each panel the posteriors for optimal estimation (red) and Bootstrap Monte 
Carlo (green) are shown. The Gaussian priors for each gas are shown with the dot dashed red curve. The true answer is the vertical black 
line. 
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Fig. 12. — Temperature profile posteriors using the Level-by-Level temperature profile approach for each observational scenario (columns) 
and tw o of the retrieval technique (rows). The solid black curve in each panel is the true temperature profile constructed with equations 
|13|16| and the parameters in Table^ The blue curve is the median temperature profile. The dashed black curve is the priortemperature 
profile constructed from Xa, as in Figure W] The prior widths for each level (not shown) are ±400 K. The dark and light red regions are 
the 1- and 2(t (68% and 95%) uncertainties, respectively. The green curve is the averaging kernel profile for temperature. The atmospheric 
regions over which this is a maximum is where we can retrieve temperature information with less dependence on the prior (see text). 
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Fig. 13. — The effect of three different Level-by-Level temperature profile priors on the retrieved temperatures (top) and spectra (bottom). 
The different temperature pro file priors are shown as the colored dashed curves. The prior widths (not shown) at each level are ±400 K, 
similar to those in Figure |12| The resultant retrieved profiles are shown as the solid colored cur ves. The thick black curve is the true 
temperature profile. The solid grey region is the la confidence interval from the retrievals in Figure [121 Note how the retrieved profiles all 
converge within the Icr confidence interval regardless of the temperature prior. The best agreement is in the mi ddle atmosphere where the 
thermal emission weighting functions are a maximum, and hence the averaging kernel profiles from Figure |12| are also a maximum. The 
spectra in the second row illustrate the effects of the different retrieved temperature profiles of corresponding color. There is virtually no 
difference in the resultant spectra for high quality data. The dotted curves at the bottom are the broadband filter functions. 



